/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __ORDINALLOGISTIC_HH__
#define __ORDINALLOGISTIC_HH__

#include <algorithm>
#include <vector>

#include "feature.hh"
#include "model.hh"
#include "normalize.hh"

namespace playerpiano
{
  class OrdinalLogisticRegression
    {
      private:
        unsigned int            n_features;
        unsigned int            n_labels;
        unsigned int            n_items;
        float*                  weights;
        std::vector<float>      logpvec;
        float                   priorw_sq_stddev;

      public:
        OrdinalLogisticRegression (unsigned int _n_features,
                                   unsigned int _n_labels,
                                   unsigned int _n_items,
                                   float*       _weights,
                                   float        priorw_stddev)
          : n_features (_n_features),
            n_labels (_n_labels),
            n_items (_n_items),
            weights (_weights),
            logpvec (n_labels),
            priorw_sq_stddev (priorw_stddev * priorw_stddev)
          {
          }

        OrdinalLogisticRegression (Model* model)
          : n_features (1U << strtoul (static_cast<char*> (model->lookup ("--n_feature_bits").first), 0, 0)),
            n_labels (strtoul (static_cast<char*> (model->lookup ("--n_labels").first), 0, 0)),
            n_items (strtoul (static_cast<char*> (model->lookup ("--n_items").first), 0, 0)),
            weights (static_cast<float*> (model->lookup ("weights").first)),
            logpvec (n_labels),
            priorw_sq_stddev (1)
          {
          }

        void
        apply_deferred_prior (unsigned int f,
                              uint64_t     s,
                              uint64_t     t,
                              float        eta,
                              float        rho)
          {
            // L2 regularization, zero mean, priorw_stddev
            // only for feature weights
            // threshold are densely updated so densely regularized
            
            if (t > s)
              {
                float powt = fastpow (t, 1.0f - rho);
                float sigma = priorw_sq_stddev;
                float pows = fastpow (s, 1.0f - rho);
                float decay = fastexp (eta * (pows - powt) / 
                                       (n_items * (1.0f - rho) * sigma));

                float* w = weights + (n_labels - 1) + (f % n_features);

                *w *= decay;
              }
          }

        // not actually ordinal logistic regression
        // polytomous rasch function form is more convenient

        const std::vector<float>&
        logp (const IndexedFeatureSet& features)
          {
            float z = 0;

            for (IndexedFeatureSetIterator it = begin (features);
                 it != end (features);
                 ++it)
              {
                unsigned int f = (*it).feature;
                float v = (*it).value;

                float* w = weights + (n_labels - 1) + (f % n_features);

                z += *w * v;
              }

            logpvec[0] = 0;

            for (unsigned int j = 1; j < n_labels; ++j)
              {
                logpvec[j] = logpvec[j - 1] + z - weights[j - 1];
              }

            // logp[j] = logpraw[j] - log (sum_i exp (logpraw[i]))
            // dlogp[j] / dlogpraw[i] 
            //  =  delta_ij - exp (logpraw[i]) / (sum_i exp (logpraw[i])) 
            //  =  delta_ij - exp (logp[i])
            //  =  delta_ij - p[i]
            // dlogpraw[0] / dw[f] = 0
            // dlogpraw[j] / dw[f] = dlogpraw[j - 1] / dw[f] + (f.v)
            //                     = j * (f.v)
            // dlogpraw[0] / dw[l] = 0
            // dlogpraw[j] / dw[l] = -1 1_{j > l}

            safe_normalize (logpvec);

            return logpvec;
          }

        // P (x) and Q (x)
        // objective function F = \sum_i q[i] logp[i]
        // so (dF / dlogp[i]) = q[i]
        // so (dF / dlogpraw[i]) 
        //   = \sum_j (dF / dlogp[j]) (dlogp[j] / dlogpraw[i])
        //   = \sum_j q[j] (delta_ij - p[i])
        //   = q[i] - p[i] (\sum_j q[j])
        // assuming q is normalized
        //   = q[i] - p[i] 
        //   
        // weird!  same result as when a linear function of the probability
        // is used with squared error ... hmmm .... (?)
        
        void
        update (uint64_t                    t,
                float                       eta,
                float                       rho,
                const IndexedFeatureSet&    features,
                const std::vector<float>&   logp,
                const std::vector<float>&   q)
          {
            float p[n_labels];
            float etat = eta * fastpow (t, -rho);

            for (unsigned int j = 1; j < n_labels; ++j)
              {
                p[j] = fastexp (logp[j]);
              }

            for (IndexedFeatureSetIterator it = begin (features);
                 it != end (features);
                 ++it)
              {
                unsigned int f = (*it).feature;
                float v = (*it).value;

                float* w = weights + (n_labels - 1) + (f % n_features);

                for (unsigned int j = 1; j < n_labels; ++j)
                  {
                    *w += etat * (q[j] - p[j]) * j * v;
                  }
              }

            float sumd = 0;

            for (unsigned int j = n_labels; j > 1; --j)
              {
                sumd -= etat * (q[j - 1] - p[j - 1]);
                weights[j - 2] -= etat * weights[j - 2] / (n_items * priorw_sq_stddev);
                weights[j - 2] += sumd;
              }
          }

        static void
        ensure (Model* model)
          {
            unsigned int n_features (1U << strtoul (static_cast<char*> (model->lookup ("--n_feature_bits").first), 0, 0));
            unsigned int n_labels (strtoul (static_cast<char*> (model->lookup ("--n_labels").first), 0, 0));
            model->ensure ("weights", 
                           ((n_labels - 1) + n_features) * sizeof (float));
          }
    };
}

#endif // __ORDINALLOGISTIC_HH__
